Operators in HySoP, basic ideas

This is an introduction about all kinds of operators available in HySoP with the basic ideas and common tools shared between all operators. You should read this whole section before any attempt to use a specific operator.

Continuous operators

A continuous operator is an object that handles a set of continuous fields, depending on time and space, with an ‘apply’ function that may modify some of these variables.

For example, consider the advection of a scalar \(\rho\) at velocity v, as defined in the following system

\[\frac{\partial{\rho(x, t)}}{\partial{t}} + v(x,t).\nabla \rho(x,t) = 0\]

To simulate this problem, an advection operator must be defined as Code

>>> from hysop import Box, Field
>>> from hysop.operators import Advection
>>> from hysop.parameters.scalar_parameter import ScalarParameter
>>> box = Box(dim=3)
>>> v = Field(domain=box, name='v',  is_vector=True)
>>> rho = Field(domain=box, name='rho',  is_vector=False)
>>> advec = Advection(v, rho,
...                   variables={v:(33, )*3, rho:(33, )*3},
...                   dt=ScalarParameter('dt'))
>>> g = advec.to_graph()
>>> _ = g.build()
>>> g.apply()

Obviously, many other parameters must be specified, such as which methods will be used for the space discretisation, the time integration, time step value and so on. This depends on each operator and will be detailed later.

An operator is determined by

  • its name which sets the ‘type’ of problem it handles,

  • a list of variables, the continuous fields, defined on one and only one domain,

  • an mpi-task,

  • some i/o parameters.

Variables

Each operator is associated with a set of continuous fields splitted into ‘input’ (read only) and ‘output’ (read and write) variables.

Computational operators

Methods

Method argument must be a dictionnary with some of the following keys/values:

Methods have default values, taken from methods preset dictionnaries.

Discretisation and setup

Each computational operator has a ‘discretize’ method, which:

and a ‘setup’ function which creates a discrete operator and ensure that the operator is ready for simulation. This function handles also internal work arrays management (see below).

So, the correct and required sequence for each computational operator is:

# build
op = Operator(...)
# discretize
op.discretize()
# make it ready for apply
op.setup(rwork=..., iwork=...)

Discrete operators

A discrete operator is the representation of the space discretisation of a continuous operator, i.e. a discrete object which applies on some discrete fields. Such an object is associated with at least one topology (but maybe several for multi-resolution operators).

Data-distribution operators

Some operators are not solving anything but help for manipulating data.

I/O operators

See :I/O management: tools and conventions in HySoP

Data communications

See MPI topologies and space discretisation.

Computational operators

The main operators are here described and explained in the context of the typical form of problem solved by the library, that is the Navier-Stokes equations in their non-dimensional velocity(\(\mathbf{u}\))-vorticity(\(\boldsymbol{\omega}\)) formulation:

\[\begin{split}\frac{\partial\boldsymbol{\omega}}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\omega} &= (\boldsymbol{\omega} \cdot \nabla) \mathbf{u} + \dfrac{1}{Re} \Delta \boldsymbol{\omega} + \nabla \times \mathbf{f}_{ext}\\ \Delta \mathbf{u} &= - \nabla \times \boldsymbol{\omega}\end{split}\]

In the above system of equations, the first one corresponds to the momentum equation with :

  • \((\mathbf{u} \cdot \nabla) \boldsymbol{\omega}\): the advection term

  • \((\boldsymbol{\omega} \cdot \nabla) \mathbf{u}\) : the stretching term (only in 3D). This term vanishes in 2D.

  • \(\dfrac{1}{Re} \Delta \boldsymbol{\omega}\) : the diffusion term with \(Re = \frac{L \mathbf{u}_{\infty}}{\nu}\), with \(L\) the characteristic length, \(\mathbf{u}_{\infty}\) the characteristic upstream/far field velocity and \(\nu\) the constant kinematic viscosity of the fluid.

  • \(\nabla \times \mathbf{f}_{ext}\) : the external forcing term that depends on the problem being solved

The second equation, \(\Delta \mathbf{u} = - \nabla \times \boldsymbol{\omega}\), is the Poisson equation allowing to recover the velocity \(\mathbf{u}\) from the vorticity \(\boldsymbol{\omega}\). This equation is derived from the incompressibility condition \(\nabla \cdot \mathbf{u} = 0\) and the definition of the vorticity field \(\boldsymbol{\omega} := \nabla \times \mathbf{u}\).

Building on this foundation, the library can handle the modeling and solving of a wide range of problems.

Numerical approach

The particularity of the HySoP library relies on the fact that the problems are solved by using a semi-Lagrangian method: the so-called “remeshed Vortex method” or “remeshed particle method”. The momentum equation can be split into transport and diffusion terms, by relying on the operator splitting methods. The idea behind the proposed numerical method is to split the equations such that each subproblem can be solved by using a dedicated solver based on the most appropriate numerical scheme and by employing a space discretization that is regular enough to be handled by accelerators (GPUs).

Semi-lagrangian (remeshed) particle methods makes it possible to solve convection problems in Lagrangian way (on particles) without imposing a Courant-Friedrichs-Lewy stability constraint, but a less restrictive stability condition allowing the use of larger time steps. In order to avoid the distortion of the convected fields, the vorticity and scalar values carried by each particle are distributed (after the advection step) on the neighboring points of an underlying Cartesian mesh. This step is called the “remeshing”.

On the one hand this method is particularly adapted for problems dominated by transport phenomena, in particular the transport of scalars in high Schmidt number flows. On the other hand, through the presence of an underlying grid, this method allows the use of eulerian solvers. In particular the library uses Cartesian grids that are compatible with a wide variety of numerical methods such as finite difference methods and spectral methods.

Operator splitting

Viscous -or operator- splitting within the framework of vortex methods was originally proposed in the early seventies by Chorin (1973). Its convergence has been proven by Beale and Majda (1982) and reformulated by Cottet and Koumoutsakos (2000) to suit to the present vorticity-velocity formulation.

The operator splitting approach is here described in the context of the Navier-Stokes equations in their velocity-vorticity formulation:

\[\begin{split}\frac{\partial\boldsymbol{\omega}}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\omega} &= (\boldsymbol{\omega} \cdot \nabla) \mathbf{u} + \dfrac{1}{Re} \Delta \boldsymbol{\omega} + \nabla \times \mathbf{f}_{ext}\\ \Delta \mathbf{u} &= - \nabla \times \boldsymbol{\omega}\end{split}\]

1. Poisson equation: PoissonCurl or Poisson

\[- \Delta \mathbf{u} = \nabla \times \boldsymbol{\omega}\]

2. Diffusion: Diffusion

\[\frac{\partial\boldsymbol{\omega}}{\partial t} = \dfrac{1}{Re} \Delta \boldsymbol{\omega}\]

3. Advection: Advection or DirectionalAdvection

\[\frac{\partial\boldsymbol{\omega}}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\omega} = 0\]

4. Stretching: DirectionalStretching

\[\frac{\partial\boldsymbol{\omega}}{\partial t} = (\boldsymbol{\omega} \cdot \nabla) \mathbf{u}\]

5. External forces: SpectralExternalForce

\[\frac{\partial\boldsymbol{\omega}}{\partial t} = \nabla \times \mathbf{f}_{ext}\]

An extensive description of the provided operators is available in API documentation of the operator package.